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I.  INTRODUCTION 


A.  BACKGROUND 

Great  research  has  been  accomplished  in  the  field  of  underwater  explosions 
(UNDEX)  since  World  War  II.  Surface  ships  and  submarines  were  subjected  to  damage 
caused  by  mines,  depth  charges  and  torpedoes  throughout  the  war.  Direct  contact  with 
such  ordinance  was  not  required  for  major  damage  to  occur.  The  requirement  to  be  able 
to  predict  and  reduce  damage  to  a  ship’s  hull  from  an  UNDEX  event  was  evident.  The 
US  Navy  has  been  performing  live  fire  ship  shock  trials  over  the  years  in  accordance  with 
OPNAVINST  9072.2  [Ref  1],  NAVSEA  0908-EP-000-3010A  [Ref  2]  and  MIE-S-901D 
[Ref  3]  to  observe  the  UNDEX  phenomenon  and  its  effects  on  existing  NAVY  ships. 

The  lead  ship  of  each  class  is  subjected  to  these  trials,  and  subsequent  trials  are 
conducted  at  the  discretion  of  the  CNO  to  validate  corrective  actions  taken  from 
deficiencies  identified  in  the  first  trial,  to  validate  changes  to  a  class  of  ship  not 
represented  in  the  first  trial,  or  to  validate  the  work  of  different  shipbuilders.  Due  to 
major  upgrades  to  the  Arleigh  Burke  guided  missile  destroyer  (DDG)  class,  the  USS 
WINSTON  CHURCHILE  (DDG-8I)  underwent  these  ship  shock  trials  in  May  and  June 
of  2001 .  Three  separate  explosions  at  various  distances  were  conducted,  and  with  the  use 
of  over  six  hundred  sensors  located  throughout  the  ship,  the  dynamic  response  to  the 
UNDEX  was  recorded. 

Computer  modeling  has  been  utilized  to  predict  many  of  the  aspects  from  an 
UNDEX  event.  Also,  the  finite  element  method  has  become  a  standard  among  engineers 
to  calculate  dynamic  responses  in  complex  structures.  Combining  these  two  theories, 
finite  element  models  of  a  ship  and  its  surrounding  fluid  have  been  subjected  to  virtual 
attacks.  These  UNDEX  simulations  have  been  remarkably  accurate  in  predicting  the 
early  time  vertical  response  when  compared  to  the  actual  ship  shock  trial  data  collected. 

With  just  a  few  keystrokes  the  geometry  of  a  charge  may  be  changed,  or  the 
physical  amount  of  charge  may  be  increased.  Without  harm  to  crew,  hull,  mechanical, 
electrical  systems  and  to  the  marine  environment  a  plethora  of  conditions  can  easily  be 
modified.  While  only  three  explosions  were  conducted  with  DDG-8 1 ,  the  uses  of  these 
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simulations  show  incredible  potential.  In  addition,  eomputer  teehnology  advances  have 
greatly  redueed  the  amount  of  time  required  to  perform  these  simulations.  However, 
there  are  still  areas  where  the  simulation  results  do  not  match  the  actual  trial  data.  In 
particular,  the  athwartship  motion  has  yet  to  be  modeled  with  as  great  success  as  the 
vertical  motion. 

In  the  near  future  these  virtual  models  will  be  more  accurate,  run  faster,  and 
perhaps  eliminate  the  need  for  live  fire  testing.  The  ability  to  predict  and  correct 
weaknesses  in  a  partieular  hull  design  before  a  ship  is  built  shows  unlimited  potential  for 
the  future  in  terms  of  cost  savings  and  operational  readiness. 


B,  SCOPE  OF  RESEARCH 

Using  the  techniques  developed  at  the  Naval  Postgraduate  School  and  with  actual 
ship  shock  trial  data  collected  in  May  and  June  of  2001  of  the  USS  WINSTON 
CHURCHILL  (DDG-8I),  this  thesis  serves  to  further  the  researeh  eonducted  by 
Sehneider  [Ref  4]  and  Didoszak  [Ref  5].  Model  refinement  and  comparative  analysis  of 
the  vertieal  and  athwartship  veloeity  response  have  been  the  basis  of  these  previous 
studies.  This  thesis  will  further  examine  the  athwartship  veloeity  response  and  present 
shock  spectra  analysis.  Data  processing  techniques  will  also  be  examined  along  with 
reasons  behind  the  disparities  that  exist  between  the  simulated  values  and  the  actual  trial 
data. 
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II.  SHIP  SHOCK  MODELING 


A.  FINITE  ELEMENT  MODEL 

Two  finite  element  models  of  Arleigh  Burke  Class  Destroyers  were  ereated  by 
Gibbs  and  Cox,  Ine.  [Ref  6]  to  assist  with  the  ship  shoek  modeling  effort.  Modifieations 
made  to  this  elass  of  ship  are  known  as  the  Flight  IIA  DDG’s.  The  USS  WINSTON 
CHURCHILL  (DDG-8I)  was  modeled  in  far  greater  detail  than  the  previous  finite 
element  model  USS  John  Paul  Jones  (DDG-53)  as  shown  below  in  Table  1. 


Model  Summary 

DDG-53 

DDG-81 

Number  of  Nodes 

35,216 

40,513 

Degrees  of  Freedom  (6  per  node) 

211,296 

243,078 

Shell  Elements 

41,078 

48,603 

Bar  and  Beam  Elements 

42,659 

49,968 

Discrete  Elements 

416 

416 

Rigid  Elements 

55 

55 

Table  1 .  Finite  Element  Model  Specifics  Comparison  [From  Ref  6.] 


Some  of  the  major  modifications  to  this  ship  class  include  the  following: 

•  Installation  of  Dual  Helicopter  hangers 

•  Extension  of  the  Transom  by  five  feet 

•  Replacement  of  the  five  inch  54  caliber  gun  with  a  five  inch  62 
cailber  gun 

•  Removal  of  forward  and  aft  VLS  cranes  and  installation  of  6  new 
VLS  cells 


•  Addition  of  five  blast  hardened  bulkheads 
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As  shown  in  Figure  1,  the  physical  differences  between  the  DDG’s  are  apparent. 
Thus,  shock  trials  needed  to  be  conducted  in  keeping  with  OPNAVISNT  9072.2  [Ref  1]. 


Figure  1 .  DDG-81  (top)  and  DDG-53  (bottom)  Finite  Element  Models  [From  Ref.  6] 

These  models  were  generated  in  a  MSC  NASTRAN  compatible  environment  and 
were  validated  by  performing  physical  checks  onboard  each  of  the  vessels.  Different 
techniques  were  used  in  the  DDG-81  model  that  served  to  improve  model  and  simulation 
accuracy,  when  compared  to  earlier  models.  For  example,  instead  of  smearing  mass 
across  the  model  the  actual  DDG-81  weight  database  (weight-tapes)  was  incorporated. 
Much  of  the  finite  element  model  was  generated  with  the  use  of  Computer  Aided 
Drawing  (CAD)  blueprints  of  the  vessel.  This  technique  enabled  the  model  to  be  built 
much  more  quickly  and  far  more  precise  than  the  DDG-53  model.  Also,  the  liquid  load, 
ordinance,  cargo  and  personnel  loads  during  each  explosion  was  calculated  and 
incorporated  into  the  model.  Since  there  were  three  explosions,  or  three  “Shots,”  three 
separate  finite  element  models  of  DDG-81  were  used  in  the  simulations. 

Figure  2,  shows  clearly  the  great  amount  of  detail  involved  with  the  DDG-81 
finite  element  model.  The  nominal  mesh  size  of  this  model  is  27  inches  by  48  inches. 
However,  Combat  Information  Center  (CIC)  was  modeled  with  a  finer  mesh  density  due 
to  the  amount  of  sensitive  electronics,  allowing  a  more  detailed  analysis  of  that  particular 
compartment. 
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Figure  2.  Centerline  eutaway  view  of  DDG-81  (port  side)  Finite  Element  Model  [From  Ref  6] 


B,  FLUID  MESH 

The  fluid  mesh  is  the  finite  element  model  of  the  water  that  surrounds  the  ship 
model.  Created  at  the  Naval  Postgraduate  Sehool  (NPS)  by  Flart  [Ref.  8]  and  Schneider 
[Ref.  4],  this  mesh  is  essential  for  the  simulation  of  an  UNDEX  event  with  LS- 
DYNA/USA  software.  Using  TrueGrid,  a  mesh-generating  program,  the  wetted  surface 
of  the  DDG-81  finite  element  model  was  extruded  into  a  finite  element  model  of  its  own. 
Figure  3  shows  the  fluid  inner  liner,  which  is  orthogonal  along  the  wetted  surface  of  the 
ship  with  the  exception  of  the  waterline  and  liner  seams.  The  liner  seams  exist  since 
areas  such  as  the  sonar  dome,  keel,  bow  and  stem  that  have  complex  shapes.  These 
seams  have  been  smoothed  by  adding  curved  wedges  to  prevent  gaps.  At  only  seven 
inches  thick  this  inner  liner  serves  as  the  foundation  for  additional  fluid  mesh  layers. 
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Figure  3 .  Thin  Fluid  Mesh  Inner  Liner  [From  Ref.  8] 


The  next  layer  of  the  fluid  mesh  liner  is  an  additional  eight  inehes  thiek.  The 
eombined  two-layer  inner  mesh  ensures  that  the  eritieal  element  thiekness  of  nine  inehes 
required  by  the  Underwater  Shoek  Analysis  Code  (USA)  is  not  exeeeded.  This  eombined 
mesh  also  establishes  the  waterline  at  21.5  feet.  Although  additional  layers  were  added  to 
this  mesh,  it  was  at  this  point  that  the  virtual  modeling  eould  begin.  The  additional  layers 
are  known  as  the  inner  mesh,  transitional  fluid  mesh,  and  outer  fluid  mesh. 

Comprehensive  studies  eondueted  by  Hart  [Ref  8]  and  Sehneider  [Ref  4]  varied 
the  ultimate  depth  of  the  outer  fluid  mesh.  Results  of  their  simulations  have  shown  that 
the  eombined  fluid  mesh  must  reaeh  the  bottom  of  the  lower  eavitation  boundary.  Figure 
4  shows  the  eoupled  model  with  the  various  layers  of  the  fluid  mesh. 
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Figure  4.  Fluid  Mesh  Layers:  Inner  Liner  (a),  Inner  Mesh  (b), 
Transitional  Mesh  (e),  Outer  Mesh  (d)  [From  Ref.  7] 


C.  NODE  IDENTIFICATION  AND  SELECTION 

Previous  simulations  in  the  Shock  and  Vibrations  Computational  Lab  (SVCL)  at 
NPS  examined  the  same  series  of  thirty-two  nodes  from  within  the  finite  element  model. 
Although  these  nodes  represented  a  variety  of  areas  of  interest  from  the  entire  ship,  this 
study  focused  more  attention  on  choosing  the  specific  nodes  that  matched  the  exact 
locations  of  the  sensors.  This  was  performed  so  that  that  the  motion  simulated  in  this 
study  could  be  compared  to  actual  trial  data  without  any  doubt. 

Having  been  provided  with  a  list  of  three-dimensional  coordinates  for  each  nodes 
in  the  DDG-81  finite  element  model  from  Gibbs  and  Cox  [Ref  10]  and  a  list  of  the  three- 
dimensional  coordinates  of  each  dynamic  sensor  from  NAVSEA  [Ref  11]  the  task  of 
matching  these  positions  took  place.  With  over  600  sensors  locations  and  over  175,000 
total  nodes  (combined  ship  and  fluid  model)  a  computer  program  was  written  to  perform 
this  task. 
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All  of  the  nodes,  sensors,  and  respeetive  coordinates  were  input  in  ASCII  format 
into  a  MATLAB  script  file.  At  this  point  the  three-dimensional  distance  formula  was 
applied  to  calculate  the  distance  between  each  sensor  with  every  single  node.  This  data 
was  stored  in  a  vector,  and  an  iterative  command  commenced  the  same  task  for  the  next 
sensor.  When  completed,  a  matrix  of  data  was  stored  holding  all  of  the  distances 
between  every  node  and  every  sensor.  MATLAB ’s  built-in  min  command  was  then  used 
to  locate  the  least  distance  between  each  sensor  and  each  node. 

With  this  data,  the  best  pair  of  nodes  and  sensors  were  matched.  There  exists 
other  methods  of  performing  this  task,  but  this  method  seemed  to  provide  the  best  results. 
For  example,  had  a  weighted  averaging  technique  been  applied  to  a  cube  of  eight  nodes 
containing  one  sensor,  the  simulations  would  have  been  adversely  prolonged.  All  of  the 
nodes  analyzed  in  this  study  were  less  than  eight  inches  away  from  the  location  of  the 
sensors.  In  comparison,  the  previous  studies  node  selection  criteria  was  based  on 
locations  within  the  nominal  mesh  size  of  a  sensor’s  location.  This  allowed  for  node 
selection  with  distances  as  great  as  48  inches  from  its  corresponding  sensor. 


D.  SIMULATING  UNDEX  WITH  LS-DYNA/USA 

LS-DYNA  is  an  industry  standard  finite  element  analysis  program  that  is  capable 
of  three-dimensional  non-linear  dynamic  analysis  using  explicit  time  integration 
algorithms.  Underwater  Shock  Analysis  Code  (USA)  is  a  series  of  complex  algorithms 
that  approximate  the  fluid  behavior  in  an  UNDEX  event.  These  two  programs  have  been 
combined  into  one  powerful  tool  that  generates  a  six  degree  of  freedom  dynamic  analysis 
to  a  floating  or  submerged  structure.  Due  to  the  fact  that  the  physics  of  underwater 
explosions  is  rather  complex,  a  concise  background  of  this  phenomenon  will  be 
described. 

1,  UNDEX  Phenomenon 

There  is  a  defined  sequence  that  occurs  during  an  UNDEX  event.  The  first  of 
which  is  the  incident  pressure  wave  or  incident  shock  wave.  This  is  known  as  a  shock 
wave  due  to  the  fact  that  the  pressure  waves  velocity  exceeds  the  acoustic  velocity  of 
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water.  The  shock  waves  energy  propagates  through  the  water  under  extreme  pressure. 
This  pressure  can  easily  be  calculated  at  any  point,  as  long  as  the  point  is  located  between 
10  and  100  times  the  charge  radius  away  from  the  charge  location.  The  USA  code 
effectively  calculates  the  incident  pressure  wave. 

When  the  shock  wave  reaches  the  free  surface,  it  becomes  reflected  as  a  tensile 
pressure  wave.  Using  an  image  charge  positioned  above  the  waterline  symmetrically 
from  the  actual  charge  the  tensile  wave  can  be  modeled.  Figure  5  shows  the  how  the 
image  charge  is  utilized.  Mathematically,  the  sum  of  the  pressure  wave  and  the  tension 
wave  would  yield  a  negative  value  over  a  large  area.  However,  since  water  does  not 
sustain  tension  cavitation  occurs.  This  is  commonly  known  as  the  bulk  cavitation  effect. 
Properly  calculating  this  bulk  cavitation  zone  is  essential  to  modeling  an  UNDEX  event. 
Figure  6  shows  the  bulk  cavitation  zone  of  a  10,000-pound  HBX-1  charge.  Local 
cavitation  also  occurs  at  the  interface  between  the  fluid  and  structure. 


DEPTH 


Figure  5.  Charge  and  Image  Charge  Positioning  [From  Ref.  9] 
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Figure  6.  Bulk  Cavitation  Zone  (not  to  scale) 


An  explosive  gas  bubble  ereated  upon  detonation  will  expand  rapidly  underwater. 
Unlike  an  explosion  in  air,  the  pressure  inside  of  the  gas  bubble  will  reaeh  equilibrium 
with  the  hydrostatie  pressure  of  the  surrounding  fluid.  The  radius  of  the  bubble  will 
naturally  inerease  until  this  equilibrium  is  reaehed,  and  momentum  will  eause  the  bubble 
radius  to  continue  to  increase  beyond  equilibrium.  At  it’s  maximum  radius  the  bubble 
will  collapse  due  to  the  positive  pressure  gradient  between  the  bubble  and  surrounding 
fluid.  At  this  point  the  bubble  is  small,  yet  a  negative  pressure  gradient  exists  between 
the  bubble  and  surrounding  fluid  causing  the  bubble  to  expand  again  towards 
equilibrium.  Depending  on  the  depth  and  size  of  the  charge  this  bubble  oscillation  will 
continue  to  occur  as  it  migrates  toward  the  surface  and  as  the  energy  is  dissipated  into  the 
surrounding  fluid. 

2.  Simulating  UNDEX  with  USA 

Using  a  Doubly  Asymptotic  Approximation  (DAA)  the  USA  code  calculates  the 
fluid  structure  interaction  during  the  UNDEX.  This  method  is  extremely  effective  at 
approximating  the  early  time  (high-frequency)  response  and  late  time  (low-frequency) 
response  [Ref.  12].  Starting  with  the  second  order  linear  equation  of  motion 

[M]{x}+[C]{x}+[K]{x}  =  {f}  (1) 
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where  [M],  [C],  and  [K]  are  the  mass,  damping,  and  stiffness  matriees 
respectively.  The  Force  vector  {f}  is  defined  as 

{f}  = -[G][Af  ]{(pj+ps)}  +  {fD}  (2) 

where  [G]  is  the  transformation  matrix  that  relates  the  structural  and  fluid  surface 
forces,  [Af]  is  a  diagonal  area  matrix  pertaining  to  the  fluid  elements,  {pi}  and  {ps}  are 
the  incident  and  pressure  wave  vectors,  and  {fo}  is  the  dry  structure  applied  force  vector. 

Having  a  defined  equation  of  motion  is  essential  to  calculating  the  dynamic 
response  to  this  or  any  finite  element  model.  DAA  utilizes  the  structural  equation  of 
motion  and  a  fluid  particle  equation  of  motion  shown  below 

[Mf]{Ps}+pc[Af]{pd  =pc[Mf]{ud  (3) 

where  [Mf]  is  the  mass  matrix  of  the  fluid  mesh,  and  {us}  is  the  scattered-wave 
fluid  particle  velocity  vector.  The  following  equation  utilizes  the  transformation  matrix 
described  in  equation  (2)  to  relate  the  fluid  particle  velocities  and  the  structural  response. 

[Gf  {x}  =  {u,}+{us}  (4) 

Combining  equations  (1)  through  (4)  the  following  two  equations  are  formed. 

[M]{x}+[C]{x}+[K]{x}  =-[G][Af]{(p,+Ps)}  (5) 

[Mf]{ps}+pc[Af]{ps}=pc[Mf]([Gf  {x}-{u,})  (6) 

The  USA  code  solves  these  systems  of  equations  for  the  forces  at  the  fluid 
structure  interface.  Knowing  these  forces,  finite  element  analysis  software  may  be  then 
used  to  calculate  the  structural  response.  Using  the  following  three  programs  within  USA 
the  shock  wave  is  generated  and  the  DAA  method  is  used  to  calculate  the  forces: 

•  FLUMAS  -  Generates  the  fluid  mass  matrix  and  the  pressure  from  the 
spherical  shock  wave  [Ref.  13] 

•  AUGMAT  -  Takes  data  from  FLUMAS  and  assembles  the  DAA  matrices 
in  preparation  for  solution  [Ref  13] 
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•  TIMINT  -  Conducts  a  step  by  step  numerieal  time  integration  using  the 
DAA  method  [Ref.  13] 

3,  Finite  Element  Analysis 

Prior  to  any  simulation  effort  the  mass,  damping,  and  stiffness  matriees  need  to  be 
determined.  The  struetural  model  and  shoek  trial  data  of  DDG-53  led  to  extensive 
researeh  eondueted  at  NPS  in  proportional  damping.  Finite  element  analysis  software 
may  be  very  powerful,  but  it  requires  linear  symmetrieal  matriees.  Unlike  the  mass  and 
stiffness  matriees,  whieh  are  easily  developed  from  material  properties,  the  damping 
matrix  must  be  approximated  so  that  the  property  of  orthogonality  holds  true.  If  it  is  not 
orthogonal  a  eonsistent  damping  eoefficient  may  not  be  used.  Using  proportional 
damping  the  matrix  [C]  from  equation  (1)  is  approximated  by  a  linear  eombination  of  the 
mass  and  stiffness  matriees  as  shown  below 

[C]  =  a[M]+p[K]  (7) 

Sinee  the  mass  and  stiffness  matriees  are  symmetric,  a  linear  eombination  of  the 
two  matriees  would  also  be  symmetrie.  The  coeffieients  alpha  and  beta  were  not 
ealeulated  with  ease,  months  of  researeh  of  actual  shock  trail  data  was  utilized  with  least 
squares  eurve  fitting  techniques  to  ealeulate  their  values.  Yet  onee  established  these 
eoeffieients  serve  to  form  an  approximation  of  the  damping  matrix  that  may  now  be 
diagonalizable  using  orthogonality.  The  eoeffieients  determined  in  DDG-53  were  used  in 
the  DDG-81  simulation  [Ref.  14]. 

LS-DYNA  was  used  to  ealeulate  the  structural  response  to  this  combined  ship  and 
fluid  model.  The  DDG-81  model  was  developed  in  a  MSC  NASTRAN  eompatible 
environment  and  subsequently  was  eonverted  into  a  keyword  format  for  LS-DYNA 
proeessing  at  NPS.  The  following  figure  depiets  the  entire  simulation  proeess. 
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Figure  7.  Simulation  Process  Flowchart  [From  Ref.  4] 


LS-DYNA  is  capable  of  performing  many  kinds  of  analysis  including  plastic 
deformation,  and  is  widely  used  in  the  automotive  industry  for  crash  test  simulations. 
However,  all  of  the  simulations  presented  in  this  study  are  using  elastic  material 
properties.  Six  degree  of  freedom  analysis  is  conducted,  and  using  post-processing  and 
viewing  software  such  as  Ceetron  Glview  nodal  displacement,  velocity,  and  acceleration 
vectors  can  easily  be  exported  and  plotted. 
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III.  ATHWARTSHIP  RESPONSE 


Although  much  smaller  in  magnitude  when  eompared  to  the  vertieal  response,  as 
shown  in  Figure  8,  the  athwartship  response  is  still  a  very  important  aspeet  to  shook  and 
vibrations  oommunity.  Muoh  of  the  highly  sensitive  eleotronio  equipment  installed 
aboard  today’s  war  fighting  ships  are  free  standing  with  isolation  mounts  only  on  the 
deok.  There  is  essentially  no  shook  proteotion  for  this  equipment  when  subjeoted  to 
athwartship  motion.  The  intent  of  this  study  is  to  validate  the  researoh  being  oonduoted  at 
NFS  in  ship  shook  modeling  by  improving  the  athwartship  motion  simulations. 


Velocity  Comparision 

Grid  222240  (A2105) 


Vertical  Response  Athwartship  Response""] 


Figure  8.  Comparison  Between  Vertieal  and  Athwartship  Veloeity  Response 

A,  DRIFT 

Measuring  the  athwartship  motion  during  the  ship  shock  trials  on  DDG-81  triaxial 
accelerometers  were  utilized.  These  sensors  capture  and  store  the  aooeleration  digitally 
in  unit’s  of  gravity.  However,  all  of  the  time  history  oomparisons  presented  are  veloeity 
vs.  time.  Due  to  that  faot,  all  of  the  sensor  data  provided  by  NAVSEA  [Refill]  needed  to 
be  integrated  and  sealed  appropriately.  Drift  oeeurring  in  the  vast  majority  of  the  vertieal 
response  integrated  data  was  reported  by  Sehneider  [Ref  4].  This  is  beeause  the  sensor  is 
sampling  a  frequency  range  greater  than  it  was  intended  for.  When  the  aeeeleration  time 
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history  data  is  integrated  the  zero  position  of  the  veloeity  eurve  shifts  away  from 
equilibrium.  To  be  certain  that  the  numerical  integration  routine  being  used  was  stable 
and  was  not  the  root  source  of  the  drift,  several  different  and  more  robust  integration 
algorithms  were  utilized,  all  with  the  same  results. 

Drift  compensation  has  been  incorporated  into  a  post  processing  software  tool, 
UERDtools.  This  routine  requires  at  least  1000  msec  of  time  history  data  to  function.  It 
was  shown  to  be  effective  for  correcting  drift  in  the  vertical  response  case  [Ref  4]. 
However,  the  athwartship  response  time  history  curves  from  the  ship  shock  trails  have 
significant  amounts  of  drift.  Using  the  drift  compensation  routine  in  UERDtools  does  not 
eliminate  the  drift.  However,  residual  drift  still  remains  in  the  athwartship  time  history 
curves  after  applying  this  technique,  leaving  the  simulation  comparisons  far  less  precise 
as  the  vertical  comparisons  [Ref.  5]. 

Developing  a  more  robust  method  of  drift  compensation  and  applying  it  to  the 

ship  shock  trial  data  was  required  to  be  able  to  conduct  an  in  depth  analysis  of  the 

athwartship  motion.  After  studying  many  of  the  ship  shock  trial  velocity  curves  it  was 

apparent  that  the  drift  resembled  second  and  third  degree  polynomial  curves.  By 

subtracting  a  second  or  third  degree  curve  from  the  time  history  plot  would  remove  the 

drift.  Using  linear  algebra  algorithm  in  MATLAB  that  solves  for  a  polynomial  curve  of 

best  fit  of  the  original  data,  and  then  calculating  the  difference  between  the  original  data 

and  the  fitted  curve  a  drift  compensated  curve  would  be  computed.  This  method  in  its 

simplest  form  described  above  is  quite  powerful,  therefore  the  user  must  observe  caution 

since  it  does  change  the  data.  Recent  upgrades  to  UERDtools  have  included  a 

polynomial  drift  compensation  function,  however  it  requires  the  specification  of  a  start 

time  in  which  the  correction  begins.  Having  to  specify  a  start  time  for  correction  requires 

judgment  that  may  leave  too  much  room  for  poor  engineering  practice.  Utilizing  the 

linear  algebra  routine  developed  in  MATLAB  requires  no  judgment  since  all  of  the  data 

is  used  in  the  curve  fitting  process.  A  compromise  can  easily  me  achieved  by  using  the 

UERDtools  polynomial  drift  compensation  function  applied  to  all  of  the  data.  This 

method  produces  the  same  results  as  the  MATLAB  function  if  and  only  if  the  time 

history  curves  begins  at  zero  msecs.  Appendix  C  includes  a  detailed  description  of  the 

linear  algebra  routine  and  a  complete  program  listing.  The  following  figure  shows  a 
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comparison  between  the  two  drift  compensation  methods  compared  with  the  original 
sensor  data. 


DDG-81  SHOX  3 


A7008A 


- A.7008A.  (Least  Squares  Method) 

- A.7008A. 

- A.7008A  (E>rift  Compensated) 


Figure  9.  Drift  Compensation  Methods 


B.  NOISE 

High  frequency  noise  tends  to  clutter  the  time  history  curves  of  the  ship  shock 
trial  data.  Frequencies  above  250  Hz  have  little  to  no  effect  to  the  overall  motion  of  the 
ship.  Applying  low  pass  filters  (filtering  all  data  above  250  Hz)  to  both  of  the  simulated 
data  and  to  the  ship  shock  trial  data  have  been  performed  at  NFS  throughout  recent 
studies  and  have  been  shown  to  greatly  improve  data  correlation  of  the  mean,  variance 
and  standard  deviations  [Ref.  4].  Figures  10  and  1 1  shows  low  pass  fdtering  applied  to  a 
set  of  simulation  data.  Note  that  the  shape  of  the  curve  is  not  changed  as  the  filter  is 
applied.  Filtering  in  this  study  has  shown  to  remove  the  small  oscillations  from  both  the 
shock  trial  data  and  the  simulation  data  allowing  for  greater  reliability  when  comparing 
curves.  Without  filtering,  there  is  a  band  of  noise  that  surrounds  the  curves  making 
statistical  comparisons  between  the  various  sensors  and  nodes  difficult  to  reproduce. 
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Grid  212058  (A.2102A.) 


I  Unfiltered  Data  | 


Figure  10.  Unfiltered  Simulation  Data 


Grid  212058  (A2102A) 


- Filtered  Data 


Figure  1 1 .  Filtered  Simulation  Data  (250  Flz  -  Low  Pass  Filter) 
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C.  ERROR  MEASUREMENT 

Russell’s  Error  Factor  [Ref.  15]  was  used  in  the  previous  DDG-81  studies  at  NFS 
and  was  very  effective  in  characterizing  the  vertical  response.  This  statistical  formula 
compares  differences  in  magnitude  and  differences  in  phase  as  two  separate  functions.  A 
comprehensive  error  factor  is  developed  from  the  square  root  of  the  sum  of  the  squares  of 
the  two  calculations.  From  this  comprehensive  error  factor,  there  are  three  categories  that 
transient  data  may  be  categorized  as:  Excellent,  Acceptable,  and  Poor.  However,  recent 
research  in  the  athwartship  motion  left  some  doubt  about  Russell’s  Error  Factor.  Since 
there  were  no  nodes  that  were  characterized  as  excellent  in  Didoszak’s  study  [Ref  5] 
perhaps  the  three  categories  do  not  accurately  depict  the  data  presented. 

The  magnitudes  of  the  athwartship  motion  are  three  to  four  times  less  than  their 
vertical  counterparts,  yet  the  simulations  are  not  as  nearly  as  precise  according  to  the 
same  statistical  criteria.  Eooking  more  closely  at  the  combined  time  history  plots  and 
evaluating  the  Russell’s  error  factor  it  can  easily  be  seen  that  the  magnitude  error  is  very 
low,  yet  as  the  time  increases  over  150  msecs  the  phase  between  the  two  curves  seems  to 
vary.  The  general  trend  of  the  curve  is  modeled  quite  well,  yet  it  appears  to  be  out  of 
phase  as  time  increases  to  500  msecs.  Figures  12  through  14  depict  this  trend  common  to 
the  majority  of  the  athwartship  simulations. 
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Figure  12.  50  msec  Comparison 


19 


1313Ci-81  SHOT  3 

C3rita  2  12058  C 1  02A) 


Figure  13.  100  msec  Comparison 
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Figure  14.  Full  500  msec  Comparison 


Russell’s  error  changes  dramatically  as  time  increases  throughout  the  simulation. 
Table  2  shows  the  relatively  low  and  consistent  magnitude  errors.  The  phase  error 
however,  which  is  extremely  low  in  the  early  time  history,  grows  beyond  the  acceptable 
range.  The  phase  error  seems  to  be  the  greatest  obstacle  to  overcome  in  the  athwartship 
simulations.  All  error  comparisons  previously  conducted  at  NFS  were  performed  with 
500  msec  time  history  curves. 
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Grid  212058  (A2102A) 

Magnitude  Error  (RM) 

Phase  Error  (RP) 

Comprehensive  Error 
(RC) 

Figure  12-50  msec 

0.15 

0.07 

0.15 

Figure  13  -  100  msec 

0.10 

0.21 

0.21 

Figure  14  -  500  msec 

0.14 

0.38 

0.36 

Table  2.  Russell’s  Error  Factor  Comparison  as  Time  Increases 
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IV.  INTRODUCTION  TO  THE  SHOCK  SPECTRA 


Previous  studies  conducted  at  NPS  have  concentrated  solely  on  time  history 
analysis  and  efforts  to  validate  the  computer  modeling  effort.  Introducing  shock  spectra 
analysis  shows  the  data  comparison  between  computer  simulations  and  shock  trials  in  the 
frequency  domain.  This  allows  for  various  aspects  of  UNDEX  and  simulations  to  be 
compared  that  aren’t  easily  recognizable  in  the  time  domain. 


A,  SHOCK  SPECTRA  BACKGROUND 

The  shock  spectra  is  defined  as  the  maximum  absolute  response  of  an  undamped 
single  degree  of  freedom  system  produced  by  a  shock  loading  [Ref  9].  If  one  were  to 
calculate  the  response  of  a  system  at  a  certain  frequency  a  curve  would  be  generated. 
Using  iterative  programming  the  response  of  a  system  can  be  characterized  by  a  series  of 
curves,  each  curve  representing  the  response  for  a  unique  frequency.  Instead  of 
analyzing  many  different  curves,  it  is  convenient  to  view  the  maximum  absolute  value 
from  each  frequency.  These  maximum  values  plotted  on  one  curve  form  the  shock 
spectra.  Time  history  plots  can  be  used  to  generate  shock  spectra  plots  with  a  simple 
algorithm.  UERDtools  has  a  very  useful  shock  spectra  generating  function  that  enables 
quick  generation  of  these  spectra  plots  in  a  variety  of  formats  including  triaxial  plots. 

Since  there  is  no  phase  in  a  shock  spectra  plot,  the  phase  concerns  described  in  the 
previous  chapter  are  eliminated.  Also,  being  able  to  pay  particular  attention  to  certain 
frequencies  shows  unlimited  potential  for  future  design  of  shock  isolation  mounts  and 
systems.  As  another  method  of  validating  the  use  of  computer  models  for  naval  ship 
shock  trials  there  seems  to  be  more  practical  uses  involved  with  the  shock  spectra 
analysis  than  with  the  time  history  analysis. 


B,  EVALUATING  THE  SHOCK  SPECTRA 

The  following  figure  is  an  example  of  a  shock  spectra  plot  and  its  corresponding 
time  history  plot  including  curves  from  shot  1  of  the  DDG-81  ship  shock  trial  and  the 
NPS  computer  simulation  of  the  same. 
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Athwartship  Velocity  (ft/s)  Pseudo  Velocity  (ft/sec) 


DDG-81  SHOX  1 

Grid  414356  (A21 13A) 


Frequency  (Hz) 

Ship  Shock:  Trial  - NPS  Simulation 


DDG-81  SHOX  1 

Grid  414356  (A21 13A) 


Xime  (msec) 

Ship  Shock  Trial  - NPS  Simulation 


Figure  15.  Sample  Shock  Spectra  Plot  (top)  and  Corresponding  Time  History  (bottom) 
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Analyzing  the  amount  of  data  presented  in  this  plot  may  be  somewhat 
overwhelming  at  first,  but  essentially  it  is  very  easy  to  understand.  First,  the  x  and  y  axes 
are  both  logarithmie.  The  y-axis  is  labeled  “Pseudo  Veloeity,”  due  to  the  fact  that  the 
peak  response  occurs  after  the  event  (UNDEX  in  this  particular  case).  Being  in  the 
frequency  domain  vice  the  time  domain  it  is  easy  to  compare  the  response  at  specific 
frequencies,  most  importantly  the  lower  natural  frequencies  of  the  structure.  The 
diagonal  and  off-diagonal  axes  also  provide  the  relative  displacement  and  acceleration. 
For  example,  to  read  the  acceleration  response  at  a  certain  frequency  we  would  identify 
the  point  at  which  the  curve  intersects  that  particular  frequency,  then  follow  the  diagonal 
axis  down  and  to  the  right  of  the  plot.  Likewise,  for  reading  the  relative  displacement 
response  at  that  same  frequency  we  would  again  start  at  the  intersection  of  the  curve  at 
that  particular  frequency,  then  follow  the  off-diagonal  axis  up  and  to  the  right  of  the  plot. 
The  top  and  right  sides  of  the  plot  include  values  for  the  displacement  and  acceleration  in 
logarithmic  form. 

Below  10  Hz  the  simulation  data  and  the  actual  data  curves  are  nearly  on  top  of 
each  other.  The  low  frequency  data  in  the  shock  spectra  plot  originates  from  the  early 
time  response.  In  this  case  the  early  time  response  of  the  simulation  is  very  accurate  as 
shown  in  corresponding  time  history  plot. 

Between  10  Hz  and  100  Hz  there  are  some  fluctuations  between  the  two  curves 
and  above  100  Hz  there  is  an  even  higher  degree  of  fluctuation.  It  is  essential  to 
recognize  that  this  is  a  logarithmic  plot  and  approx  two-thirds  of  all  of  the  data  presented 
is  below  1  ft/sec.  For  example,  there  appears  to  be  a  gap  between  the  simulation  and 
actual  data  at  100  Hz,  however  there  is  only  0.3  ft/sec  difference  between  the  two  curves. 
When  compared  to  the  peak  response  of  the  vertical  motion,  0.3  ft/sec  may  be  considered 
negligible.  Without  using  the  log  scale,  the  two  curves  would  essentially  be  blurred 
together  and  reading  information  from  the  plots  would  be  difficult.  Yet,  the  log  scale 
does  prove  without  any  doubt  that  the  simulation  effectively  models  the  low  frequency 
response  of  the  UNDEX  attack  extremely  well. 
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V.  EVALUATION  AND  ANALYSIS 


A.  TIME  HISTORY  ANALYSIS 

Each  time  history  from  the  ship  shoek  trial  athwartship  data  used  in  the  DDG-8 1 
eomparison  had  signifieant  drift.  As  deseribed  in  Chapter  III,  polynomial  eurve  fitting 
was  utilized  to  remove  this  drift.  In  previous  studies  Russell’s  Error  Eaetor  has  been  used 
to  benehmark  the  level  of  aeeuraey  in  the  NFS  simulations.  An  analysis  eondueted  an  of 
the  athwartship  response  and  eoneluded  that  the  drift  assoeiated  with  the  sensor  data  from 
the  aetual  ship  shoek  trial  overwhelmed  the  rather  small  magnitudes  eomputed  in  the 
simulation  [Ref  5.].  Also  noted  in  the  study  was  that  the  late  time  phase  error  in 
eombination  with  the  signifieant  drift  resulted  in  poor  data  eorrelation  using  the  Russell’s 
error  criteria. 

Eigure  17,  is  an  athwartship  response  eomparison  of  Russell’s  Error  Eaetor  from  a 
previous  study  [Ref  5].  All  of  the  data  eompiled  in  this  figure  was  taken  from  500  msec 
time  history  plots.  Eigure  18  shows  the  Russell’s  Error  Eaetor  eomparisons  from  shot  1 
in  this  study.  The  pink  squares  depiet  the  Russell’s  Error,  both  phase  and  magnitude  for 
the  100  msec  time  history.  It  is  clear  from  the  figures  that  there  is  a  marked  improvement 
in  study  presented  in  this  paper  due  to  the  reduetion  of  the  phase  error.  More  than  half  of 
points  have  a  magnitude  error  of  less  than  0.1  that  ean  most  likely  be  attributed  to  the 
polynomial  drift  eompensation  teehnique. 

The  vast  majority  of  the  phase  error  is  from  the  last  half  of  the  time  history. 
Although,  the  magnitude  of  the  response  may  be  aecurately  simulated  the  eurves  are 
often  eompletely  out  of  phase  after  the  first  250  msees.  It  is  in  this  late  time  portion  of 
time  history  that  the  response  is  typieally  less  than  0.2  ft/see.  The  peak  response  in  is 
eaptured  in  the  first  70-80  msees  in  the  three  shots  analyzed,  and  the  simulations  have 
performed  exceptionally  well  in  modeling  the  early  time  response.  Yet,  aecording  to  the 
Russell’s  Error  Eaetor  guidelines  the  phase  and  magnitudes  errors  are  weighed  equally  in 
determining  the  eomprehensive  error.  Because  it  seems  that  the  negligible  late  time 
athwartship  motion  is  adversely  affecting  the  integrity  of  the  simulation  efforts  a  similar 
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comparison  was  performed  using  just  the  first  100  msecs.  The  eomparison  made  with  the 
100  msee  time  history  is  much  better  than  the  500  msee  eomparison. 


Similar  Results  were  notieed  in  the  analysis  of  shot  2  as  shown  in  Figure  19  and 
Shot  3  as  shown  in  Figure  20.  Shot  2  has  historieally  been  the  least  aeeurate  athwartship 
simulation  based  on  Russell’s  Error  Faetor  comparisons  [Ref.  5].  However,  when 
evaluating  just  the  first  100  msees  the  data  eorrelates  well  with  two  eomparisons  in  the 
exeellent  range  and  the  majority  in  the  aeceptable  range. 


Russell's  Comprehensive  Error  Factor 
Athwartship  Analysis  (Shots  1,  2  &  3) 
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Figure  16.  Russell’s  Error  Factor  Comparison  from  Previous  Study  [From  Ref  5] 


28 


Phase  Error  Phase  Error 


DDG-81  SHOT  1 

Russell's  Error  Factor  Comparison  of  Athwartship  Velocities 


Figure  17.  Russell’s  Error  Factor  Comparison  for  Shot  1 


DDG-81  SHOT  2 

Russeii's  Error  Factor  Comparison  of  Athwartship  Veiocities 


Figure  18.  Russell’s  Error  Factor  Comparison  for  Shot  2 
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DDG-81  SHOT  3 

Russell's  Error  Factor  Comparison  of  Athwartship  Velocities 


Figure  19.  Russell’s  Error  Factor  Comparison  for  Shot  3 


Outliers  within  the  shock  trial  data  have  been  previously  identified  at  NFS  [Ref. 
4].  Having  evaluated  a  different  set  of  sensors  and  nodes  in  this  study  specific  bow  and 
stern  sensors  identified  in  the  past  were  avoided.  However,  when  comparing  the  nodal 
response  from  the  simulations  there  were  a  few  instances  in  which  the  actual  sensor  data 
seemed  corrupt.  For  example,  Figure  20  shows  one  of  these  outliers.  In  every  case 
evaluate  the  peak  response  occurred  within  the  first  100  msecs.  However,  this  sensor 
data  has  captured  a  peak  response  of  -6  ft/sec  at  approximately  350  msecs.  The 
amplitude  of  the  peak  response  of  this  sensor  raises  doubts  as  to  its  accuracy.  Compared 
with  hundreds  of  other  athwartship  response  plots,  this  sensor’s  response  does  not  seem 
to  dampen.  This  sensor  when  compared  to  simulation  results  had  a  high  Russell’s  Error 
Factor  in  all  three  shots. 
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DDG-81  SHOX  3 

Grid  500048  (A2232A) 


Ship  Shock  Trial  NPS  Simulation  | 


Figure  20.  Example  of  Outlier 


B.  SHOCK  SPECTRA  ANALYSIS 

As  described  earlier  there  is  no  phase  in  a  shock  spectra  plot.  Therefore  the 
concerns  with  late  time  phasing  between  the  simulation  and  the  actual  trial  data  have 
been  eliminated.  Evaluating  the  data  in  the  frequency  domain  allows  for  a  different 
perspective  on  the  physical  behavior  of  the  UNDEX  attack  and  the  simulation.  As  seen 
from  various  time  history  plots,  the  magnitudes  of  the  athwartship  motion  are  small  in 
comparison  to  the  magnitudes  of  the  vertical  motion.  This  fact  is  also  depicted  in  the 
shock  spectra  plots.  The  majority  of  all  the  data  presented  in  the  athwartship  shock 
spectra  plots  is  below  1  ft/sec  in  velocity.  Some  of  the  peak  values  shown  in  this  study 
are  between  2  and  3  ft/sec  and  can  be  compared  to  a  gentle  swell  in  the  ocean. 

NPS  simulations  have  performed  exceptionally  well  in  the  1  to  10  Hz  range,  there 
is  little  difference  between  the  two  curves  in  this  range.  The  majority  of  these  plots 
shown  in  Appendix  B.  display  a  gradual  rise  in  amplitude  as  the  frequency  increases. 
Between  10  and  50  Hz  the  peak  response  has  occurred  and  there  is  a  downward  trend. 
However,  there  tends  to  be  some  spikes  between  10  and  20  Hz.  Above  50  Hz  the 
response  fluctuates  a  great  deal  but  the  downward  trend  is  prevailing.  It  can  be  shown 
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that  the  NPS  simulations  over  prediet  the  high  frequeney  responses  from  approximately 
70  Hz  upwards.  These  values  all  fall  below  1  ft/see,  and  on  the  log  seale  the  eurves  still 
remain  elose.  The  upper  limit  for  all  of  these  shoek  speetra  plots  was  set  at  250  Hz  since 
the  data  was  low-pass  filtered.  The  following  table  is  a  general  summary  of  the  Shock 
Spectra  analysis  by  frequency  range. 


Frequency 

Trend 

Shot  1 

Shot  2 

Shot  3 

1  to  10  Hz 

Gradual  rise  to 
peak  value 
(approx.  3  ft/see) 

Simulation  elosely 
matehes  trial  data 

Simulation  elosely 
matehes  trial  data 

Simulation 
routinely  under 
prediets  trial  data 

10  to  30  Hz 

Oseillation  near 
peak  values 
(within  1  to  2 
ft/see) 

Simulation  elosely 
matehes  or  slightly 
over  prediets  trial 
data 

Simulation  elosely 
matehes  or  slightly 
over  prediets  trial 
data 

Simulation  leads 
trials  data  by  5  Hz 

30  to  70  Hz 

Consistent 
response  below 
peak 

Simulation  elosely 
matehes  or  slightly 
over  prediets  trial 
data 

Simulation 
routinely  over 
prediets  trial  data 

Simulation 
routinely  over 
prediets  trial  data 

70  to  100  Hz 

Peak  value  oeeurs 
(approx.  3  ft/see) 
then  sharp  deeline 

Simulation 
routinely  over 
prediets  trial  data 

Simulation 
routinely  over 
prediets  trial  data 

Simulation 
routinely  over 
prediets  trial  data 

100  to  250  Hz 

High  degree  of 
oseillation  (below 

1  ft/see) 

Simulation  elosely 
matehes  or  slightly 
over  prediets  trial 
data 

Simulation  elosely 
matehes  or  slightly 
over  prediets  trial 
data 

Simulation  elosely 
matehes  or  slightly 
over  prediets  trial 
data 

Table  3 .  Summary  of  Shock  Spectra  Analysis  of  Shots  1 ,  2  and  3 


C.  NUMERICAL  EVALUATION  OF  DATA 

Throughout  the  course  of  this  study  great  attention  has  been  paid  to  concepts  such 
as  drift  and  late  time  phase  disparities.  The  sensors  used  in  the  DDG-81  ship  shock  trial 
were  12-bit  sensors  calibrated  to  record  data  between  -200g  and  200g.  This  range  has 
adequately  recorded  the  vertical  motion  yet  12  bit  technology  handles  2  (4096)  possible 

values.  Therefore,  the  step  size  is  fixed  with  a  resolution  0.098g  which  is  calculated  by 
dividing  the  calibration  span  by  the  step  size.  In  other  words,  the  sensor  is  limited  in  its 
ability  to  accurately  capture  the  high  frequency,  low  amplitude  response  that  is  typical  of 
the  athwartship  motion  [Ref.  16].  Reducing  the  span  that  the  sensors  are  calibrated  to  a 
lower  value  would  increase  the  resolution  of  the  late  time  response.  This  would  also  be 
very  impractical  since  the  vertical  response  would  essentially  be  chopped. 
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For  example,  if  16  bit  sensors  were  utilized  instead  of  12  bit  sensors  during  the 
ship  shoek  trial  the  amount  steps  available  for  reeording  would  be  2'^  (65536).  Therefore 
the  resolution  would  be  0.006 Ig  using  the  same  calibration  span.  Both  the  vertical  and 
athwartship  responses  could  be  accurately  recorded  and  the  amount  of  drift  shown  in  the 
late  time  response  would  be  reduced  or  even  eliminated.  However,  high  frequency  noise 
filtering  may  become  even  more  important  when  processing  data  since  the  ability  to 
capture  this  range  is  improved  by  an  order  of  magnitude. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  SUMMARY  OF  FINDINGS 

Although  very  low  in  magnitude  when  compared  to  the  vertical  response,  the 
athwartship  response  remains  an  important  aspect  of  an  UNDEX  attack.  Being  able  to 
accurately  predict  the  athwartship  response  serves  as  a  means  validate  this  technique. 
This  thesis  has  focused  on  the  reasons  why  past  simulations  have  yet  to  model  the 
athwartship  response  with  the  degree  of  accuracy  as  the  vertical  response  and  ways  to 
improve  the  simulation  results. 

The  largest  obstacle  in  modeling  the  athwartship  response  is  the  reliability  of 
actual  ship  shock  trial  data.  12  bit  sensors  used  in  the  trials  have  not  adequately  recorded 
the  late  time  response  leaving  significant  amounts  of  drift  in  the  integrated  velocity  data. 
Reducing  this  drift  has  been  the  key  to  improving  the  accuracy  of  the  simulations. 
Russell’s  Error  Eactor  is  a  technique  used  to  compare  two  sets  of  data.  However,  if  a  set 
of  actual  trial  data  is  corrupt  the  simulation  data  will  correspond  poorly.  In  other  words, 
the  simulation  may  be  accurate,  but  the  error  will  be  high.  Eooking  beyond  statistical 
error  is  essential  when  evaluating  the  differences  between  the  simulation  and  trial  data. 

Using  more  strict  criteria  in  the  selection  of  nodes  has  been  shown  to  modestly 
improve  the  simulation  results.  These  improvements  have  not  been  nearly  as  substantial 
as  those  shown  by  compensating  the  drift  with  polynomial  curve  fitting.  However, 
selecting  the  node  that  best  represents  the  location  of  a  sensor  is  always  desirable. 

Shock  spectra  analysis  can  serve  as  a  design  tool  as  well  as  a  tool  for  comparative 
analysis.  Barge  testing  has  been  used  to  shock  qualify  naval  equipment  for  years,  yet 
using  these  UNDEX  simulations  and  the  shock  spectra’s  created  accurate  predictions  of 
the  frequency  response  can  be  achieved.  As  a  comparative  tool,  the  shock  spectra  has 
shown  that  the  low  frequency  response  is  very  accurately  modeled,  and  in  many  cases  the 
simulation  are  more  conservative  than  the  actual  trial  data. 

Determining  the  amount  of  data  required  to  perform  the  analysis  of  the 
athwartship  motion  is  essential.  Unlike  the  vertical  response,  500  msec  time  histories  of 
the  athwartship  response  may  be  too  much  information.  When  evaluating  just  the  first 
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100  msecs  the  peak  response  is  captured  and  has  been  simulated  to  the  level  of  accuracy 
of  previous  vertical  response  studies.  The  late  time  response  is  so  low  in  magnitude  that 
it  may  be  considered  negligible.  Gentle  sea  swells  greatly  exceed  the  average  late  time 
response.  Clearly,  a  response  that  low  in  magnitude  is  not  going  to  be  a  concern  of  any 
surface  ship  commander,  naval  architect,  marine  engineer,  or  anyone  else  involved  the 
UNDEX  community. 


B,  FUTURE  WORK 

Being  able  to  analyze  different  types  of  UNDEX  attacks  are  of  great  interest  to 
the  military.  This  thesis  along  with  previous  studies  conducted  at  NFS  has  served  to 
benchmark  the  validity  of  using  simulations  instead  of  live  fire  shock  testing.  In  light  of 
recent  events  including  the  USS  Cole  suicide  attack  (2000),  and  the  USS  Samuel  B. 
Roberts  mine  explosion  (1988)  the  ability  to  simulate  hull  damage  due  to  close  in 
explosions  may  have  the  potential  to  save  lives.  The  whipping  effects  due  to  the  gas 
bubble  oscillation  have  been  blamed  for  sinking  many  liberty  ships  in  WWII.  This  is 
another  area  that  modeling  may  serve  to  improve  hull  design. 

Using  the  same  techniques  described  in  this  thesis,  the  UNDEX  simulation  of  a 
ship  while  still  in  the  design  phase  may  be  the  next  step.  More  specifically,  the  next 
generation  destroyer  (DDX)  can  be  modeled  in  the  same  fashion  that  DDG-81  was 
modeled.  Environmental  regulations  are  only  going  to  become  more  strict,  and  the  day 
may  come  that  live  fire  ship  shock  testing  will  be  outlawed.  These  simulations  may  serve 
as  the  only  way  to  evaluate  the  response  of  a  surface  ship  to  and  UNDEX  attack. 

In  the  near  future  the  application  of  16  bit  accelerometers  as  compared  to  the  12 
bit  accelerometers  that  were  used  in  the  DDG-81  ship  shock  trials  will  increase  the 
resolution  of  the  ship’s  response.  Increased  resolution  will  lead  to  the  elimination  of  the 
velocity  drift,  enabling  for  more  accurate  comparisons  of  these  simulations  to  the  trial 
data. 
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APPENDIX  A.  DDG-81  ATHWARTSHIP  TIME  HISTORY  PLOTS 
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Figure  2 1 .  Combat  Information  Center  (1-126-0-C)  Deek  Sensor  FM  126 
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Figure  22.  Combat  Information  Center  Annex  (1-126-0-C)  Deck  Sensor  FM  174 
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Figure  23.  Radar  Room  #1  (03-128-0-C)  Bulkhead  Sensor 
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Figure  24.  Passage  Way  (02-133-1-L)  Overhead  Sensor  FM  142 


HHO-81  SHOT  1 


Grid  4143  70  C A2  1  1  5 A.) 


Time  (msec) 

I  —  Ship  Shook.  Trial  T*«JF*S  Siini-il£itic>i~i  | 


Figure  25.  Radar  Room  #2  (03-142-0-C)  Deck  Sensor 
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Figure  26.  Radar  Room  #2  (03-142-0-C)  Deck  Sensor  FM  174 
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Figure  28.  Stbd  Mast  Leg  (MAST)  Forward  Outboard  Corner 
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Figure  29.  Central  Control  Station  (1-268-0-C)  Deck  Sensor 
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Figure  30.  Combat  Systems  Equipment  Room  #1  (2-053-1-C)  Deck  Sensor 
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Figure  3 1 .  Combat  Systems  Maintenance  Center  (01-130-0-Q)  Deck  Sensor 
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Figure  32.  Combat  Information  Center  (1-126-0-C)  Deek  Sensor  FM  126 
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Figure  33 .  Combat  Information  Center  Annex  (1-126-0-C)  Deck  Sensor  FM  174 
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Figure  34.  Radar  Room  #1  (03-128-0-C)  Bulkhead  Sensor 
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Figure  35.  Passage  Way  (02-133-1-L)  Overhead  Sensor  FM  142 
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Figure  36.  Radar  Room  #2  (03-142-0-C)  Deck  Sensor 
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Figure  37.  Radar  Room  #2  (03-142-0-C)  Deck  Sensor  FM  174 
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Figure  38.  Port  Mast  Leg  (MAST)  Forward  Outboard  Corner 
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Figure  39.  Stbd  Mast  Leg  (MAST)  Forward  Outboard  Corner 
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Figure  40.  Central  Control  Station  (1-268-0-C)  Deck  Sensor 
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Figure  41 .  Combat  Systems  Equipment  Room  #1  (2-053-1-C)  Deck  Sensor 
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Figure  42.  Combat  Systems  Maintenance  Center  (01-130-0-Q)  Deck  Sensor 
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Figure  44. 


Combat  Information  Center  Annex  (1-126-0-C)  Deck  Sensor  FM  174 
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Figure  45.  Radar  Room  #1  (03-128-0-C)  Bulkhead  Sensor 
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Figure  46.  Passage  Way  (02-133-1-L)  Overhead  Sensor  FM  142 
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Figure  47.  Radar  Room  #2  (03-142-0-C)  Deck  Sensor 
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Figure  48.  Radar  Room  #2  (03-142-0-C)  Deck  Sensor  FM  174 
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Figure  49.  Port  Mast  Leg  (MAST)  Forward  Outboard  Corner 
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Figure  50.  Stbd  Mast  Leg  (MAST)  Forward  Outboard  Corner 
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Figure  5 1 .  Central  Control  Station  (1-268-0-C)  Deck  Sensor 
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Figure  52.  Combat  Systems  Equipment  Room  #1  (2-053-1-C)  Deck  Sensor 
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Figure  53.  Combat  Systems  Maintenance  Center  (01-130-0-Q)  Deck  Sensor 
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Figure  55. 
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Figure  56.  Radar  Room  #1  (03-128-0-C)  Bulkhead  Sensor 
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Figure  57.  Passage  Way  (02-133-1-L)  Overhead  Sensor  FM  142 
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Figure  58.  Radar  Room  #2  (03-142-0-C)  Deck  Sensor 
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Figure  59.  Radar  Room  #2  (03-142-0-C)  Deck  Sensor  FM  174 
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HHO-81  SHOT  1 


Cirici  41^2^9  (A2241  A) 


Frequency  (Hz) 

I  —  Ship  Shook.  T~i~ial  NFS  Siini-ilatioii  | 


Figure  6 1 .  Stbd  Mast  Leg  (MAST)  Forward  Outboard  Corner 
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Figure  62.  Central  Control  Station  (1-268-0-C)  Deck  Sensor 
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Figure  63.  Combat  Systems  Equipment  Room  #1  (2-053-1-C)  Deek  Sensor 


SHOT  1 

Cirici  2  1  29  1  V  (^^42  1  8A) 


5) 


Frequency  (H!z) 

-  Ship  Shook  Trial  NTS  Simi-ilal 


Figure  64.  Combat  Systems  Maintenance  Center  (01-130-0-Q)  Deck  Sensor 
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Figure  65.  Combat  Information  Center  (1-126-0-C)  Deek  Sensor  FM  126 
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Figure  66.  Combat  Infonnation  Center  Annex  (1-126-0-C)  Deek  Sensor  FM  174 
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Figure  67.  Radar  Room  #1  (03-128-0-C)  Bulkhead  Sensor 
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Figure  68.  Passage  Way  (02-133-1-L)  Overhead  Sensor  FM  142 
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Figure  69.  Radar  Room  #2  (03-142-0-C)  Deck  Sensor 
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Figure  70.  Radar  Room  #2  (03-142-0-C)  Deck  Sensor  FM  174 
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Figure  72.  Stbd  Mast  Leg  (MAST)  Forward  Outboard  Corner 
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Figure  73.  Central  Control  Station  (1-268-0-C)  Deck  Sensor 
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Figure  74.  Combat  Systems  Equipment  Room  #1  (2-053-1-C)  Deck  Sensor 


HHO-81  SHOT  2 


.0 


Grid  2  1  29  1  V  C^42  1  8A.) 


Figure  75.  Combat  Systems  Maintenance  Center  (01-130-0-Q)  Deck  Sensor 
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Figure  77. 


Combat  Information  Center  Annex  (1-126-0-C)  Deek  Sensor  FM  174 
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Figure  78.  Radar  Room  #1  (03-128-0-C)  Bulkhead  Sensor 
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Figure  79.  Passage  Way  (02-133-1-L)  Overhead  Sensor  FM  142 
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Figure  80.  Radar  Room  #2  (03-142-0-C)  Deck  Sensor 
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Figure  8 1 .  Radar  Room  #2  (03-142-0-C)  Deck  Sensor  FM  174 
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Figure  82.  Port  Mast  Leg  (MAST)  Forward  Outboard  Corner 
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Figure  83 .  Stbd  Mast  Leg  (MAST)  Forward  Outboard  Corner 
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Figure  84.  Central  Control  Station  (1-268-0-C)  Deck  Sensor 
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Figure  85.  Combat  Systems  Equipment  Room  #1  (2-053-1-C)  Deck  Sensor 
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Figure  86.  Combat  Systems  Maintenance  Center  (01-130-0-Q)  Deck  Sensor 
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APPENDIX  C.  MATLAB  DATA  PROCESSING  ROUTINES 


A.  NODE-SENSOR  IDENTIFICATION 

The  following  program  code  was  written  using  MATLAB.  The  purpose  of  this 
program  is  to  compute  the  distance  between  nodes  and  sensors,  then  to  select  a  node¬ 
sensor  pair  that  is  closest  to  each  other. 

%LT  Doug  Petrusa,  USCG 

%This  program  loads  into  memory  the  x,y,  and  z  position  of  each  node  and 
%each  sensor.  It  then  calculates  the  distance  formula  for  between  each 
%node  and  each  sensor  and  stores  that  information  in  a  Matrix  "d."  Using 
%the  "min  function"  in  MATLAB  the  smallest  distance  between  the  two  is 
%selected  and  stored  in  "locate"  and  "match".  Post  Processing  in  Excel  is 
%needed  since  the  Sensor  numbers  need  to  be  added, 
clear  all;close  all;clc; 
load  sensor  node; 
tic 

%  Node=rand(2500,4);  smaller  matrix  to  trouble  shoot  code 
%  Sensor=rand(  100,4);  smaller  matrix  to  trouble  shoot  code 
d=zeros(length(Node),length(Sensor)); 
for  i=l  :length(Node) 
for  j=  1 :  length(Sensor) 

d(iJ)=sqrt((Node(i,2)-Sensor(j,2))'^2+(Node(i,3)-Sensor(j,3))^2+(Node(i,4)-Sensor(j,4))^2); 

end 

end 

[closest  index]=min(d); 

locate(:,2)=closesf; 

locate(:,l)=index’; 

A=locate(:,l); 

B=Node; 
for  p=l:length(A) 
match=B(A,:); 
end 

match(:,5)=locate(:,2); 

toe 
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B,  POLYNOMIAL  DRIFT  COMPENSATION 

In  an  effort  to  improve  data  eorrelation  between  eomputer  simulation  and  aetual 
experimental  results  we  have  been  using  the  Drift  Compensation  Funetion  built  into 
UERD  Tools.  This  method  has  proven  effeetive  to  remove  signifieant  drift  from  the 
sensors  time  history  plot. 

However  in  many  eases  the  Drift  Compensation  funetion  in  UERDtools  seems  to 
drag  the  final  data  point  in  the  time  history  plot  baek  to  equilibrium.  Also,  this  funetion 
requires  a  minimum  of  1000  msee’s  to  operate. 

Using  simple  Einear  Algebra  and  Statistieal  Theory  I  have  written  a  very  simple 
algorithm  whieh  implements  least  squares  polynomial  eurve  fitting.  This  fitted  eurve 
represents  the  aetual  drift  or  error.  Then  subtraeting  the  drift  eurve  away  from  the  sensor 
a  eorreeted  time  history  eurve  is  obtained. 

In  Mathematieal  terms: 

The  basie  Einear  Algebra  Equation: 

Ax  =  b 

Where  A  is  a  matrix  and  x  and  b  are  eolumn  veetors  respeetively.  Substituting 
this  equation  in  terms  of  x  and  y  data  from  time  history  plots,  b  beeomes  the  y  (response) 
veetor  and  x  remains  as  the  x  (time-independent)  veetor.  Matrix  A  will  serve  as  a 
funetion  operator,  therefore  our  new  equation  is  Ax  =  y.  In  most  applieations  x  is  solved 
by  a  number  of  different  algorithms,  (eg  A'^b,  EU,  QR,  rref. . .). 

In  least  squares  eurve  fitting  we  do  not  solve  the  above  equation.  Instead  we 
modify  the  terms  and  generate  an  A  matrix  from  our  x  data.  In  bloek  matrix  terms  A  = 
the  first  eolumn  veetor  of  the  A  matrix  is  the  individual  terms  of  the  original 
X  veetor  squared  (a*^^^  =  x.^).  The  seeond  eolumn  veetor  of  A  is  the  original  x  veetor  (a*^^^ 
=  x).  The  third  eolumn  veetor  of  the  A  matrix  is  a  veetor  of  ones  (a*^^^  =  i).  The  b  veetor 
still  remains  unehanged. 
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Since  the  Matrix  has  three  columns  and  the  solution  vector  is  a  column  of 
coefficients  for  a  degree  polynomial  our  solution  vector  will  have  three  unknowns. 
Our  polynomial  has  the  standard  coefficients  of  p(x)  =  ax  +  bx  +  c.  Therefore  let’s  name 
our  new  unknown  vector  p.  The  newly  formed  equation  is  now: 

Ap  =  y 

Unfortunately  this  equation  can’t  be  solved  on  most,  if  not  all  occasions  due  to 
singularity.  Therefore  we  will  use  the  normal  equations  of  Linear  Algebra  to  get  a  least 
squares  solution.  This  equation  looks  like  the  following: 

A'^Ap  =  A^y 

By  left  multiplying  each  side  of  the  equation  with  the  transpose  of  the  A  matrix 
we  will  now  be  left  with  a  much  smaller  equation.  Dimensionally  this  equation  now 
takes  the  shape  of  a  (3  by  3)  (3  by  1)  =  (3  by  1)  which  is  very  easy  to  solve  and  leaves  us 
with  the  desired  polynomial  eoefficients. 

With  the  new  polynomial  coefficients  evaluated  at  every  time  step  (evaluated 
along  the  x  vector)  a  trend  will  be  displayed.  By  subtraeting  this  trend  from  the  original 
time  history  response  vector  a  drift  compensated  time  history  veetor  is  generated. 

This  algorithm  ean  easily  be  converted  to  fit  a  3’^‘*,  4**^,  5**' . ..n*’’  degree  polynomial 
by  adding  columns  to  the  A  matrix.  Also,  the  other  advantage  is  that  only  three  data 
points  are  needed,  not  1000.  The  following  program  code  was  written  in  MATLAB. 
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%LT  Doug  Petrusa 
%NPS  Shock  Team 

%Drift  Compensation  using  Linear  Algebra  Least  Squares  Polynomial 
%Curve  Fitting. 

sensor_data=textread('VA7008A.eu.asc’);  %Sensor  A7008A  integrated  and  converted  to  ASCII 
drift_data=textread('UERD_data.asc');  %Data  Drift  Compensated  in  UERD  Tools 
sensor_data2=textread('V2009VLasc');%Sensor  V2009VI  converted  to  ASCII 
drift_data2=textread('V2009VI_drift.asc');  %Data  Drift  Compensated  in  UERD  Tools 
%Break  Data  in  x  anmd  y  components  -  not  necessary  but  allows  for  easy  understanding  of  code 

x=sensor_data(:,l);  %  x  vector  -  Time  in  msec 
y=sensor_data(:,2);  %  y  vector  -  Velocity  response  in  (ft/sec) 

%Create  the  Matrix  needed  for  Ax=b,  Standard  Linear  Algebra  Equation 
A=[x.^2  X  ones(length(x),l)]; 

%Program  could  easily  be  made  a  cubic  fitting  program  by  adding  this  line 
%of  code  instead  of  the  above  line: 

%A=[x.^3  x.^2  X  ones(length(x),l)]; 

%Solve  for  polynomial  coeff  s  using  the  Normal  Equations  (A'*A)x=A'b 
%Coeff  s  to  fit  the  Standard  Polynomial  equation:  p(x)=  ax'^2  +  bx  +  c 
coeff=(  A'  *  A)\(  A'  *y) ; 

%Evaluate  the  newly  created  polynomial  at  each  data  point 
curve_fit=polyval(coeff,x); 

%Subtract  away  the  curve 

drift_comp=(y-curve_fit);  %This  is  the  newly  compensated  curve 

%Plotting  Commands 
figure(l) 

plot(x,drift_comp,'r') 
grid  on; 

xlabel(Time  (msec)') 
ylabel(’Velcity  (ft/sec)') 
title('Drift  Compensation  Methods') 
hold  on 
plot(x,y,'k') 

plot(drift_data(:,l),dnft_data(:,2),'b') 

legend('Least  Squares  Method', 'Sensor  Data',  'UERD  Tools  Algorithm', 0) 
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APPENDIX  D.  UERD  TOOLS  PROCESSING  COMMANDS 


This  appendix  contains  a  brief  description  of  several  of  the  commands  and 
functions  used  to  expedite  the  post  processing  of  the  time  history  plots  and  the 
corresponding  shock  spectra  plots.  UERDtools  is  a  powerful  post  processing  software 
tool  created  by  Mr.  Paul  Mantz  at  NSWC  Carderock,  MD.  This  program  is  a  collection 
of  popular  algorithms  frequently  used  by  mechanical  engineers  (e.g.  numerical 
integration,  unit  conversion,  unit  scaling,  fast  fourier  transforms)  in  a  Window’s  based 
environment.  Without  the  use  of  this  software,  each  of  the  steps  taken  to  post  process  the 
data  collected  would  have  required  a  separate  Fortran  or  MATLAB  code.  The 
Input/Output  interface  to  this  program  enabled  fast  and  easy  access  to  all  of  these 
algorithms. 


A,  COMMANDS 

The  following  is  a  description  of  the  major  command  needed  to  perform  data  post 
processing. 

•  Calculus  -  This  command  was  utilized  to  integrate  accelerometer  data 
from  the  ship  shock  trial  to  velocity  data  by  means  of  the  trapezoid 
method.  Differentiation  is  also  accomplished  by  finite  differencing. 

•  Scale  -  This  command  is  extremely  helpful  in  converting  between  units 
such  as  inches  and  feet  quickly  in  one  step.  Conversions  from  English  to 
SI  units  are  also  accomplished  in  the  same  manner.  Also,  the  dependant 
axis  can  be  scaled  from  seconds  to  milliseconds. 

•  Combine  Graphs  -  When  generating  each  of  the  time  history  curves 
presented  in  this  thesis  there  were  initially  two  separate  plots.  The  first 
plot  was  typically  the  Ship  Shock  Trial  data,  and  the  second  plot  was  the 
NFS  Simulation  data.  This  command  would  combine  the  two  curves  onto 
a  single  plot  in  one  step.  Automatically  the  plot  was  rescaled  to  properly 
display  both  curves,  a  legend  is  displayed  at  the  bottom  with  distinct  curve 

colors  that  can  be  preset  by  the  user. 
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•  Filter  -  Data  filtering  by  means  of  two-pole  Bessel  functions  was  also 
used  in  the  generation  of  the  time  history  plots.  In  each  plot  shown  in 
Appendix  A,  a  low  pass  filter  was  set  at  250Hz.  In  other  words,  all 
frequencies  above  250  Hz  were  eliminated  from  the  time  history.  Due  to 
this  fact  no  data  is  shown  above  250  Hz  in  the  shock  spectra  plots  in 
Appendix  B.  This  technique  as  described  in  Chapter  III,  removed  the  high 
frequency  noise  from  the  individual  curves. 

B,  BATCH  FILE  FUNCTIONS 

Most  of  the  commands  available  in  the  UERDtools  program  can  be  automated  by 
the  use  of  a  batch  file  operation.  Instead  of  performing  tasks  such  as  the  generation  of 
one  shock  spectra  plot  from  one  time  history  plot,  a  batch  file  operation  can  create 
multiple  plots  simultaneously.  By  choosing  a  batch  file  command  and  selecting  multiple 
time  history  plot  files  in  a  computer  directory,  the  same  command  will  be  performed  to 
each  file. 

The  following  batch  file  operations  files  were  used  to  create  all  of  the  plots 
presented  in  this  thesis: 

•  Calculus  -  Convert  Ship  Shock  Trial  data  from  Acceleration  (g’s)  to 
Velocity  (ft/sec). 

•  Drift  Compensation  -  Corrects  drift  in  Ship  Shock  Trial  data  (1000  msec 
time  history  required).  Note:  a  different  drift  compensation  method 
described  in  Chapter  III  was  used  in  the  preparation  of  the  plots  presented 
in  this  thesis. 

•  Linear  Interpolate  -  Convert  Simulation  output  (LS-DYNA)  from  ASCII 
text  to  UERDtools  graph  file  with  equal  time  steps  (4e-6  seconds). 

•  Scale  -  Appropriately  scale  Simulation  data  x-axis  from  secs  to  msecs  and 
y-axis  from  inches  to  feet. 

•  Shift  Time  -  Used  to  match  the  start  times  of  the  Shock  Trial  Data  curve 
and  Simulation  curve  (typically  80  msec  shift). 
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•  Sync  Start  Time  -  Used  after  eurves  are  eombined  on  to  one  plot  sinee  the 
time  steps  vary  between  the  two  eurves. 

•  Trim  Data  -  Cut  time  histories  from  2000  msees  to  500  msees. 

•  Filtering  -  Low  Pass  250  Hz  eurves 

•  Find  Russell’s  Error  -  Creates  a  text  files  with  the  Russell’s  Error  Faetor 
and  Gear’s  Error  of  eaeh  time  history  plot. 

•  Shock  Spectra  -  Generates  shoek  speetra  plots  from  time  history 
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